rm(list = ls())
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      warning = FALSE,
                      fig.align = 'center',
                      fig.width = 8,
                      fig.height = 6,
                      dev = 'jpeg',
                      dpi = 300)
#XQuartz is a mess, put this in your onload to default to cairo instead
options(bitmapType = "cairo") 
# (https://github.com/tidyverse/ggplot2/issues/2655)
# Lo mapas se hacen mas rapido
library(here)
library(kableExtra)
library(ggthemes)
library(ggrepel)
library(ggridges)
#analisis
library(ggpubr)
library(easystats) # multiples unciones analiticas
library(sf)
library(tidyverse, quietly = TRUE)
library(modelsummary)
library(terra) # replace raster
library(TropFishR)
library(mixR)
library(CCAMLRGIS)

Background

Methodology

A important piece of information for a stock evaluation refers to the biological components such as average sizes and weights across areas and years. To do this, we will explore the biological data and prepare the output to add it into stock assessment integrate model (Methot & Wetzel, 2013).

The object ohbio2 come from data exploration analysis in data request CCAMLR data. This objetc have bio information from krill.

#cargo objeto
meta <- get("METADATA")
c1 <- get("C1")
ohbio <- get("OBS_HAUL_BIOLOGY")

Join data set with master as c1 set. This join is trought obs_haul_id variable to get geoposition variables

ohbio2 <- left_join(c1, ohbio, by="obs_haul_id")
names(ohbio2)
##  [1] "c1_id.x"                   "obs_haul_id"              
##  [3] "obs_logbook_id.x"          "obs_haul_number"          
##  [5] "haul_number.x"             "vessel_name"              
##  [7] "vessel_nationality_code"   "fishing_purpose_code"     
##  [9] "season_ccamlr"             "target_species"           
## [11] "asd_code"                  "trawl_technique"          
## [13] "catchperiod_code"          "date_catchperiod_start"   
## [15] "datetime_set_start"        "datetime_set_end"         
## [17] "datetime_haul_start"       "datetime_haul_end"        
## [19] "datetime_timezone"         "depth_gear_set_end_m"     
## [21] "depth_gear_haul_start_m"   "depth_bottom_set_end_m"   
## [23] "depth_bottom_haul_start_m" "latitude_set_end"         
## [25] "longitude_set_end"         "latitude_haul_start"      
## [27] "longitude_haul_start"      "gear_type_code"           
## [29] "gear_type"                 "mesh_code"                
## [31] "trawl_net_number"          "notes"                    
## [33] "trawl_duration_depth_h"    "trawl_duration_total_h"   
## [35] "krill_greenweight_kg"      "c1_id.y"                  
## [37] "obs_logbook_id.y"          "haul_number.y"            
## [39] "taxon_code"                "taxon_scientific_name"    
## [41] "taxon_family"              "maturity_stage"           
## [43] "sex_code"                  "length_total_cm"          
## [45] "greenweight_kg"

Firsts glance. Test how many register have by year. In this case, length_total_cm by season ccamlr. Same exercise in date period date_catchperiod_start to separate dates.

ohbio3 <- ohbio2 %>%
  mutate(Year = year(date_catchperiod_start),
         Month = month(date_catchperiod_start),
         Day = day(date_catchperiod_start)) %>% 
  #toupper() para convertir los valores a mayúsculas
  mutate(sex_code = toupper(sex_code))

Save data further analysis

length481 <-ohbio3 %>% 
  dplyr::select(7, 9, 11, 12, 14, 24, 25, 29, 42, 44, 46, 47, 43, 45) |>   filter(asd_code=="481")
#save(length481, file = "length481.RData")

First thing is get different rater layer to join krill data length according different porpoises.

# Cargo linea de costa
coast <- load_Coastline()
coast1<- st_as_sf(coast) 
coast2 = st_transform(coast1, "+proj=latlong +ellps=WGS84")
# Uso las agrupaciones de Strata
strata <- st_read("~/DOCAS/Mapas/Antarctic_SHPfiles/Strata.shp",
                quiet=T)
strata=st_transform(strata, "+proj=latlong +ellps=WGS84")

Test Strata SSMU, just to know another way to join data, but this kind of spatial structuration is deprecated to mamagemente use (Figure@ref(fig:maptest2).

Grouping bio data into stratas

names(length481)
##  [1] "vessel_nationality_code" "season_ccamlr"          
##  [3] "asd_code"                "trawl_technique"        
##  [5] "date_catchperiod_start"  "latitude_set_end"       
##  [7] "longitude_set_end"       "gear_type"              
##  [9] "maturity_stage"          "length_total_cm"        
## [11] "Year"                    "Month"                  
## [13] "sex_code"                "greenweight_kg"
ohbio6 <- st_as_sf(length481 %>% 
                     drop_na(latitude_set_end), 
                   coords = c("longitude_set_end", "latitude_set_end"),  
                  crs = "+proj=latlong +ellps=WGS84")

Study Area

# y testeo el mapa
ssmap <- ggplot()+
  geom_sf(data = strata |> 
            filter(ID != "Outer") |>  
           mutate(ID = str_replace_all(ID, "Extra", "GS")),color="red")+
  geom_sf(data = ohbio6 |>
            drop_na(length_total_cm) |> 
            filter(Year>2017), 
          aes(color = length_total_cm)) +
  geom_text_repel(data = strata |>  
            filter(ID != "Outer") |>  
           mutate(ID = str_replace_all(ID, "Extra", "GS")), 
            aes(x = Labx, y = Laby, 
                label = ID), 
            min.segment.length = 0,
                             box.padding = 2,
                             max.overlaps = 10)+
  geom_sf(data = coast2, colour="black", fill=NA)+
  scale_color_viridis_c(option = "G",
                        name="Length (cm)")+
  ylim(230000, 2220000)+
  xlim(-3095349 , -1858911)+
  # coord_sf(crs = 32610)+ #sistema de prpyecccion para campos completos
  coord_sf(crs = 6932)+
  theme_bw()
ssmap
Strata Maps in 48.1

Strata Maps in 48.1

Length composition by Strata CCAMLR to visualization first. First step is group data into to poligons strata.

strata <- st_make_valid(strata)
sf4 <- st_join(strata, ohbio6)
# Save an object to a file
saveRDS(sf4, file = "sf4.rds")

Load RData

# Restore the object
sf4 <- readRDS("sf4.rds")

Data

histogram length data to viz in anthor way.

jzstrata <- ggplot(sf4 %>% 
                     mutate(ID = if_else(ID == "Extra", "GERLASHE", ID)) %>% 
               filter(Year>2010,
                      ID !="Outer"),
             aes(x=length_total_cm, 
                 y = as.factor(Month), 
                 fill=ID))+
  geom_density_ridges(stat = "binline", bins = 30, 
                      scale = 1.9, 
                      draw_baseline = FALSE,
                      alpha=0.9)+
  facet_grid(Year~ID) +   
  geom_vline(xintercept = 3.6, color = "red")+
  scale_x_continuous(breaks = seq(from = 1, to = 10, 
                                  by = 2))+
  scale_y_discrete(breaks = seq(from = 1, 
                                to = 12, by = 4))+
  scale_fill_viridis_d(name="Strata",
                       option="F")+
  theme_few()+
  xlab("Length (cm.)")+
  ylab("")
jzstrata

Now, must filter the DF

sf5 <- sf4 |> 
  dplyr::select(c("ID",
        "date_catchperiod_start",
        "length_total_cm",
        "Year",
        "Month",
        "sex_code",
         "greenweight_kg")) |> 
  mutate(ID = if_else(ID == "Extra", "GERLASHE", ID)) |>  
  filter(ID !="Outer") |> 
  data.frame()

Calculating the proportion of records per year and per stratum in a table. This table allows us to establish criteria for defining the strata on which we will analyze krill parameter. We leave out just JOIN.

# Calcular la proporción
proporcion <- (round(table(sf5$ID) / sum(table(sf5$ID))*100,2))

# Crear una tabla con las proporciones
tabla_proporcion <- as.data.frame(proporcion)
tabla_proporcion$Categoria <- rownames(tabla_proporcion)
rownames(tabla_proporcion) <- NULL

# Renombrar las columnas
colnames(tabla_proporcion) <- c("Strata", "%")

# Mostrar la tabla utilizando kbl()
kbl(tabla_proporcion, 
    caption = "Proportion length register by Strata")  |> 
  kable_classic(full_width = F, 
                html_font = "Cambria") |> 
  kable_styling(bootstrap_options = "striped", 
                latex_options = "striped")
Proportion length register by Strata
Strata % NA
BS 62.71 1
EI 6.58 2
GERLASHE 12.13 3
JOIN 0.14 4
SSIW 18.43 5

Parameters estimation k and L_inf_

Method ELEFAN ((Tobias2017?))

sf5fil <- sf5 |>
  filter(ID %in% c("BS", "EI", "GERLASHE", "SSIW")) |>  
  mutate(date_catchperiod_start = as.Date(date_catchperiod_start)) |>  
  mutate(yearly_Group = floor_date(date_catchperiod_start, "year")) |>   # New column
  drop_na(length_total_cm) 

# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLASHE", "EI")
# Crear una lista para almacenar los resultados de cada ID
lfq_results <- list()
# Iterar sobre cada ID
for (id in unique(sf5fil$ID)) {
  # Obtener el nombre correspondiente
  nombre <- nombres[which(unique(sf5fil$ID) == id)]
  # Filtrar los datos para el ID actual
  df_id <- sf5fil %>% filter(ID == id)
  # Crear un objeto lfq para el ID actual
  lfq <- lfqCreate(data = df_id,
                   Lname = "length_total_cm",
                   Dname = "yearly_Group",
                   bin_size = 0.1)
  # Agregar el resultado a la lista
  lfq_results[[id]] <- lfq
  # Graficar el objeto lfq
  plot(lfq, Fname = "catch",
       main = nombre)
}

Identified distribution

ei <- Bhattacharya(lfq_results$EI)
## Interactive session needed for Bhattacharya.
ss <- Bhattacharya(lfq_results$SSIW)
## Interactive session needed for Bhattacharya.
bs <- Bhattacharya(lfq_results$BS)
## Interactive session needed for Bhattacharya.
gs <- Bhattacharya(lfq_results$GERLASHE)
## Interactive session needed for Bhattacharya.

ahora damos asignaciones a curvas para representación

# Crear una lista para almacenar los resultados de PW_results <- list()

# Definir los nombres correspondientes a cada objeto lfq_results
nombres <- c("SSI", "BS", "GERLASHE", "EI")

# Iterar sobre cada objeto lfq almacenado en lfq_results
for (i in seq_along(lfq_results)) {
  lfq <- lfq_results[[i]]
  
  # Obtener el nombre correspondiente
  nombre <- nombres[i]
  
  # Restructurar el objeto lfq
  lfqbin <- lfqRestructure(lfq, MA = 3, addl.sqrt = TRUE)
  
  # Graficar el objeto lfq reestructurado
  plot(lfqbin, hist.col = c("white", "black"),
       image.col = c(rep(rgb(1,0.8,0.8),1000), "white", 
                     rep(rgb(0.8,0.8,1),1000)),
       ylim = c(0,max(lfqbin$midLengths+0.5)),
       main = nombre)
  
  # Ajustar curvas al objeto lfq reestructurado
  tmp <- lfqFitCurves(lfqbin, par = list(Linf=6.5, 
                                          K=0.45,
                                          t_anchor=0.5),
                      draw = TRUE, col=4, lty=2)
}

Proceed to fit Modal Progression analysis with ELEFAN

Parameters to SSIW

# Response surface analyss
res_RSAsswi <- ELEFAN(lfq = lfq_results$SSIW,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$SSIW$dates[3],
                  cross.midLength = lfq_results$SSIW$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAsswi <- ELEFAN_SA(lfq_results$SSIW, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.011048 maxTime = 30
## Emini is: -0.2661270789
## xmini are:
## 5.328667012 0.8705196602 0.2853670418 0.4288634807 0.1844906621 
## Totally it used 30.011068 secs
## No. of function call is: 1092
# run ELEFAN with genetic algorithm
res_GAsswi <- ELEFAN_GA(lfq_results$SSIW, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

Table with different method

# show results
RSAsswi <- unlist(res_RSAsswi$par)
GAsswi <- unlist(res_GAsswi$par)
SAsswi <- unlist(res_SAsswi$par)
parsswi <- round(rbind(RSAsswi, GAsswi, SAsswi),3)

parsswi %>%
  kbl(booktabs = T,
      position="ht!",
    caption = "Parametres LH to SSWI") %>%
  kable_paper("hover", 
              full_width = F)%>%
  kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE,
                font_size=9)
Parametres LH to SSWI
Linf K t_anchor C ts phiL
RSAsswi 5.448 0.640 0.790 0.500 0.500 1.279
GAsswi 5.879 0.596 0.331 0.463 0.475 1.314
SAsswi 5.329 0.871 0.285 0.429 0.184 1.393

Parametrs to EI

# Response surface analyss
res_RSAei <- ELEFAN(lfq = lfq_results$EI,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$EI$dates[3],
                  cross.midLength = lfq_results$EI$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAei <- ELEFAN_SA(lfq_results$EI, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.035548 maxTime = 30
## Emini is: -0.2111109108
## xmini are:
## 5.863710068 0.9813904597 0.873453021 0.9993372839 0.5819024704 
## Totally it used 30.035574 secs
## No. of function call is: 1329
# run ELEFAN with genetic algorithm
res_GAei <- ELEFAN_GA(lfq_results$EI, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

# show results
RSAei <- unlist(res_RSAei$par)
GAei <- unlist(res_GAei$par)
SAei <- unlist(res_SAei$par)
parei <- round(rbind(RSAei, GAei, SAei),3)
parei %>%
  kbl(booktabs = T,
      position="ht!",
    caption = "Parametres LH to Elephand Island") %>%
  kable_paper("hover", 
              full_width = F)%>%
  kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE,
                font_size=9)
Parametres LH to Elephand Island
Linf K t_anchor C ts phiL
RSAei 5.879 0.470 0.750 0.500 0.500 1.211
GAei 5.737 0.595 0.960 0.618 0.528 1.292
SAei 5.864 0.981 0.873 0.999 0.582 1.528

Parametrs to BS

# Response surface analyss
res_RSAbs <- ELEFAN(lfq = lfq_results$BS,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$BS$dates[3],
                  cross.midLength = lfq_results$BS$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAbs <- ELEFAN_SA(lfq_results$BS, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.019744 maxTime = 30
## Emini is: -0.3184663901
## xmini are:
## 5.310221989 0.8441457965 0.2786562045 0.3598986657 0.9663543715 
## Totally it used 30.019766 secs
## No. of function call is: 1190
# run ELEFAN with genetic algorithm
res_GAbs <- ELEFAN_GA(lfq_results$BS, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

# show results
RSAbs <- unlist(res_RSAbs$par)
GAbs <- unlist(res_GAbs$par)
SAbs <- unlist(res_SAbs$par)
parbs <- round(rbind(RSAbs, GAbs, SAbs),3)

parbs %>%
  kbl(booktabs = T,
      position="ht!",
    caption = "Parametres LH to BRANSFIELD STRAIT") %>%
  kable_paper("hover", 
              full_width = F)%>%
  kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE,
                font_size=9)
Parametres LH to BRANSFIELD STRAIT
Linf K t_anchor C ts phiL
RSAbs 5.448 0.640 0.790 0.500 0.500 1.279
GAbs 5.297 0.879 0.428 0.654 0.516 1.392
SAbs 5.310 0.844 0.279 0.360 0.966 1.377

Parametrs to GERLASHE

# Response surface analyss
res_RSAgs <- ELEFAN(lfq = lfq_results$GERLASHE,  
                  MA = 5,
                  Linf_range = seq(4.5, 7, 
                                   length.out = 30),
                  K_range = exp(seq(log(0.1),
                                    log(2),
                                    length.out = 30)),
                  C = 0.5,
                  ts = 0.5,
                  method = "cross",
                  cross.date = lfq_results$GERLASHE$dates[3],
                  cross.midLength = lfq_results$GERLASHE$midLengths[5],
                  contour = TRUE,
                  add.values = FALSE,
                  agemax = 5,
                  hide.progressbar = TRUE)
## Optimisation procuedure of ELEFAN is running. 
## This will take some time. 
## The process bar will inform you about the process of the calculations.
# run ELEFAN with simulated annealing
res_SAgs <- ELEFAN_SA(lfq_results$GERLASHE, 
                    SA_time = 60*0.5, 
                    MA = 5, 
                    agemax = 5,
                    seasonalised = TRUE, addl.sqrt = TRUE,
                    init_par = list(Linf = 4.5, 
                                    K = 0.5, 
                                    t_anchor = 0.5,
                                    C=0.5, 
                                    ts = 0.5),
                    low_par = list(Linf = 4,
                                   K = 0.01,
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5,
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1))
## Simulated annealing is running. 
## This will take approximately 0.5 minutes.
## timeSpan = 30.003085 maxTime = 30
## Emini is: -0.290322014
## xmini are:
## 5.624068536 0.9826401418 0.7312699817 0.3292821422 0.6789067537 
## Totally it used 30.003104 secs
## No. of function call is: 1593
# run ELEFAN with genetic algorithm
res_GAgs <- ELEFAN_GA(lfq_results$GERLASHE, 
                    MA = 5, 
                    seasonalised = TRUE, 
                    maxiter = 10, 
                    agemax = 5,
                    addl.sqrt = TRUE,
                    low_par = list(Linf = 4.5, 
                                   K = 0.01, 
                                   t_anchor = 0, 
                                   C = 0, 
                                   ts = 0),
                    up_par = list(Linf = 6.5, 
                                  K = 1, 
                                  t_anchor = 1,
                                  C = 1, ts = 1),
                    monitor = FALSE)
## Genetic algorithm is running. This might take some time.

# show results
RSAgs <- unlist(res_RSAgs$par)
GAgs <- unlist(res_GAgs$par)
SAgs <- unlist(res_SAgs$par)
pargs <- round(rbind(RSAgs, GAgs, SAgs),3)

pargs %>%
  kbl(booktabs = T,
      position="ht!",
    caption = "Parametres LH to GERLASHE") %>%
  kable_paper("hover", 
              full_width = F)%>%
  kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE,
                font_size=9)
Parametres LH to GERLASHE
Linf K t_anchor C ts phiL
RSAgs 5.448 0.640 0.790 0.500 0.500 1.279
GAgs 5.405 0.746 0.470 0.575 0.499 1.338
SAgs 5.624 0.983 0.731 0.329 0.679 1.492

Method NO LINEAR MODEL

We use firts a descomposition methot with library mixR and then we use nlmer (Bates et al. (2015)) to calcularte parametrs

# Agregar la nueva columna con los intervalos de longitud total
sf6 <- sf5 |> 
  drop_na(length_total_cm) |> 
  filter(ID %in% "GERLASHE",
         Year==2018,
         Month==5)

sf5a <-  sf5 |> 
  drop_na(length_total_cm ) |> 
  filter(ID %in% "BS",
         Year==2018)
  
# compruebo el n de comp

s_weibull = select(sf5a$length_total_cm, ncomp = 2:6, family = 'weibull')

s_normal = select(sf5a$length_total_cm, ncomp = 2:6)




# fit a Normal mixture model
mod1 = mixfit(sf5a$length_total_cm, ncomp = 4,
              pi = c(0.5, 0.5, 0.5, 0.5), 
              mu = c( 3, 4, 5, 6), 
              sd = c(0.2, 0.2, 0.2, 0.2))
mod2_weibull = mixfit(sf6$length_total_cm, family = 'weibull', ncomp = 3)

s_weibull = select(sf5a$length_total_cm, ncomp = 2:6, family = 'weibull')

s_normal = select(sf5a$length_total_cm, ncomp = 2:6)
plot(s_weibull)
plot(s_normal)

# plot the fitted model# plot the fitted model# plot the fitted model
plot(mod1)
plot(mod2_weibull)

# fit a Normal mixture model (equal variance)
mod1_ev = mixfit(sf6$length_total_cm, ncomp = 2, ev = TRUE)

References

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
Methot, R. D., & Wetzel, C. R. (2013). Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research, 142, 86–99. https://doi.org/10.1016/j.fishres.2012.10.012